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Abstract 

We study the possibility of improving the constraints on the lifetime of gravitino 
dark matter in scenarios with bilinear ii-parity violation by estimating the amount 
of cosmic-ray antideuterons that can be produced in gravitino decays. Taking into 
account all different sources of theoretical uncertainties, we find that the margin 
of improvement beyond the limits already set by cosmic-ray antiproton data are 
quite narrow and unachievable for the next generation of experiments. However, 
we also identify more promising energy ranges for future experiments. 
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1 Introduction 


Most current cosmic-ray detectors have the ability to determine the mass of the incoming 
particles - and in some cases also their charge. However, some isotopes in the cosmic- 
ray particle spectrum still remain to be observed - among those antideuterons. Donato, 
Fornengo and Salati were the hrst who suggested that antideuterons would be an inter¬ 
esting species for the search of dark matter (DM) via cosmic rays [1]. In the past 16 years 
since this pioneering work, various authors have followed up on this idea [2-13]. Indeed, 
antimatter atoms like antideuterium cannot exist in stars. Therefore, antideuterons - the 
corresponding nuclei - should not exist as astrophysical primary cosmic rays.^ The only 
astrophysical process that is expected to produce it is the spallation of cosmic rays off 
the interstellar medium (ISM). These so-called secondary antideuterons, for kinematic 
reasons, cannot be produced at kinetic energies lower than few GeV per nucleon [14]. 
If, however, DM annihilations or decays lead to the production of antideuterons, their 
spectrum would not exhibit this kinematic threshold. This leads to the conclusion that 
the observed antideuteron flux from DM annihilation or decay in the Galactic halo could 
be orders of magnitude higher than the astrophysical background at very low energies [1]. 

This theoretical argument motivated the construction of dedicated instruments 
for low-energy antideuteron searches like the General AntiParticle Spectrometer 
(GAPS) [15], which after a successful balloon-flight of a prototype in mid-2012 
(pGAPS [16,17]) is expected to fly again in the future with its hnal design. Gurrent 
multi-purpose cosmic-ray experiments, like AMS-02 and BESS, are also looking for this 
type of particles at somewhat higher energies. 

Even though it has been understood that because of tertiary production [18]^ and 
energy losses taking place during the cosmic-ray propagation [3] the difference between 
secondary and primary fluxes may not be as high as originally expected, antideuterons 
are still considered one of the most interesting species for DM searches. More recently, 
this interest has been extended to other yet unobserved cosmic-ray species like antihe¬ 
lium nuclei from DM annihilations or decays [19,20]. The growing activity of the held 
expresses itself through the organisation of dedicated events such as the recent antideu¬ 
teron workshop held at the University of Galifornia, Los Angeles. A summary paper on 
the status of the held has recently been completed [21]. 

In the present work, we investigate the case of gravitino DM within the framework of 
bilinear i?-parity violation [22,23].^ In this type of scenario the gravitino would be long- 
lived enough to be the DM in the universe but would eventually decay and produce - 
among other species - antideuterons. A main motivation for this scenario is that it leads 

^Primary cosmic rays are particles that are accelerated in the ISM, for instance by a supernova 
remnant. 

^Tertiary production refers to non-annihilating inelastic interactions of cosmic-ray antideuterons with 
the ISM leading to a migration of antideuterons from the high-energy part to the low-energy part of the 
spectrum. 

^The gravitino could also be a viable DM candidate in scenarios with trilinear i?-parity violation [24, 
25]. Antideuteron signals in this theoretical framework were studied in [12,13]. 
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to a consistent cosmological scenario, explaining the baryon asymmetry of the nniverse 
via thermal leptogenesis and avoiding any cosmological gravitino problems [23]. For a 
more detailed introduction of the model, please consult our previous work, where we 
studied constraints on the gravitino lifetime from cosmic-ray antiprotons in the same 
theoretical framework [26]. 

The structure of this paper is as follows: In the next section, we briefly revise deuteron 
and antideuteron formation in the coalescence model, derive the coalescence momentum 
from collider data, and simulate antideuteron spectra for the gravitino decay channels. In 
section 3, after studying the cross sections at stake, we compute the propagated cosmic- 
ray antideuteron fluxes at Earth, taking into account various sources of uncertainties. In 
section 4, we discuss the prospects for the detection of gravitino decays via antideuterons, 
before we come to our conclusions. In an appendix to this paper, we present some cross 
section parametrisations relevant for antideuteron propagation. 


2 Deuteron and Antideuteron Formation 

2.1 The Coalescence Model 

Deuterons and antideuterons may form in particle physics processes through final state 
interactions in a nucleon shower [27]. A usual approach to describe their formation is 
the phenomenological coalescence model [28,29]. In this prescription, independent of the 
details of the microscopic formation mechanism, a deuteron is formed when a proton and 
a neutron come sufficiently close in momentum space, i.e. when the absolute value of the 
difference of their four-momenta is below a threshold coalescence momentum: 

\Pp-Pn\<P0- (2.1) 

In order to determine the value of the coalescence momentum relevant for antideuteron 
production in gravitino decays we will make use of the antideuteron production rate 
measured by the ALEPH experiment at the LEP collider [30] . Although there are several 
measurements of antideuteron production in collider experiments - which in general do 
not lead to a common value for po [10,21] - we restrict to the ALEPH data since the 
gravitino decay channel '03/2 Zu produces antideuterons via Z boson fragmentation 
exactly as in the ALEPH measurement. Also the channels 03/2 —)■ Wi and 03/2 —)• hu are 
expected to be more similar to the case of Z boson fragmentation than to the production 
of antideuterons in pp collisions [31-33], e~p collisions [34], or collisions on or nearby 
the Upsilon resonance [35,36]. 

A common approximation for deuteron formation is that the distributions of neutrons 
and protons in momentum space are spherically symmetric and uncorrelated [27]. This 
leads to an energy distribution of deuterons that can be directly calculated from the 
product of the individual spectra of neutrons and protons [4,6] : 

dNg _ ^ TUd _ 1 _ 

dTd 6 rUpTrin ^Tj + 2mdTd dTp dT^ ’ 
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where Tp = Tn = Td/2 are the kinetic energies of protons, neutrons and deuterons. 
Therefore, this approximation is called factorised coalescence. 

As pointed out in [6], this approximation is qualitatively wrong and significantly 
underestimates the deuteron yield in high-energetic processes. This is due to the fact that 
the distributions of neutrons and protons are actually neither spherically symmetric nor 
uncorrelated. For instance, in the decay of a DM particle into a Z boson and a neutrino 
the probability for the formation of a deuteron in the fragmentation of the Z boson should 
be independent of the DM mass. This is due to the fact that the Z boson fragmentation 
process is always the same as viewed from the Z boson rest frame. However, the factorised 
coalescence approximation gives a lower yield of deuterons for larger DM masses, since 
the protons and neutrons are distributed over a larger phase space for higher injection 
energies. Another qualitatively wrong behaviour of this approximation is the possibility 
that protons and neutrons from distinct DM decays form a deuteron. This is due to the 
fact that the spectra of protons and neutrons are simply multiplied, while in principle the 
coalescence condition on the four-momenta of protons and neutrons should be applied 
on an event-by-event basis. 

This can be achieved in a Monte Carlo simulation of the decay process [6]. For in¬ 
stance, using an event generator like Pythia 6.4 [37] one can simulate the hadronisation 
of massive gauge and Higgs bosons and then apply, event by event, the coalescence 
condition on the protons and neutrons.'^ This method leads to plausible results, e.g. 
the deuteron yield in the DM decay to final states including W, Z or Higgs bosons is 
independent of the DM mass. However, this strategy also requires a lot of computing 
time to generate smooth spectra as only one deuteron or antideuteron is produced in 
(P(IO^) fragmentation processes. Moreover, the authors of [38] point out that baryon pair 
distributions are not used to tune event generators and therefore their predictions for 
deuteron production are not without uncertainty. In fact, the authors of [9] hnd a dis¬ 
crepancy of up to a factor 2-4 when comparing deuteron production between the Monte 
Carlo generators Herwig-1--1- [39,40] and Pythia 8 [41]. 

Due to the deuteron’s binding energy of 2.224 MeV [42], the coalescence of free pro¬ 
tons and neutrons is forbidden since energy and momentum cannot be conserved at the 
same time if no further particles are involved in the process. However, in a hadronic 
shower with many particles, it is easily conceivable that the overall process conserves 
energy and momentum and thus allows for deuteron coalescence.® When simulating deu- 

'^In fact, there are some combinatoric ambiguities in the deuteron coalescence in Pythia. Namely, 
it is possible that a proton fulfils the coalescence criteria with two different neutrons or vice versa. We 
checked numerically that the number of events where this is the case is very low. The probability increases 
with po, but even for po = 250 MeV less than 0.1% of the deuteron events are affected. Therefore, we 
conclude this not to be a significant effect. In our analysis we always combine the first pair of protons 
and neutrons fulfilling the coalescence criterion according to the particle order in the event record. 

^While low-energy deuteron formation proceeds via the process pn —>■ dy, this is not necessarily 
the case in high-energy nuclear collisions. In their pioneering work on deuteron formation, Butler and 
Pearson considered an explicit interaction with the optical potential of the target nucleus to take care 
of energy and momentum conservation [27]. However, this model was later found to be in conflict with 
experimental data [43,44]. The factorised coalescence model does not touch the issue at all since it does 


4 



teron coalescence in a Monte Carlo event generator, the fonr-momenta of protons and 
nentrons are given explicitly and we have no way of realistically treating the mnlti-body 
process leading to denteron formation. The most common approach in previous works 
on antideuteron signals from DM annihilation or decay is thus to consider momentum 
conservation and to determine the denteron kinetic energy from the proton and neutron 
momenta and the denteron mass of 1.8756 GeV [42]: 



(2.3) 


We will also follow this approach since it turns out that the uncertainty on the spectra 
introduced by this ambiguity is marginal for our discussion.® 

More recently, it was realised that the condition \pp — p„| < po is not not enough to 
guarantee a physically meaningful coalescence prescription in Monte Carlo event gener¬ 
ators [10]. Some of the antiprotons and antineutrons emerging in DM annihilations or 
decays are generated by the decay of metastable mother particles like the baryons A®, 
A“, A®, and Dj]', and the mesons and 5®, which have lifetimes of 

(!1(10“^®-10“^^) s [37,42]. By contrast, the nuclear interactions leading to antideuteron 
formation take place on femtometer length scales, corresponding to lifetimes of 0(10“^^) s 
for relativistic particles. Therefore, it is necessary to include a condition on the spatial 
separation of protons and neutrons as well. 

To accommodate this, Ibarra and Wild excluded weakly decaying baryons with life¬ 
times r > Imm/c from the decay chain [10].' Although this condition removes the 
most relevant metastable mother particles, namely A® and several other meta¬ 
stable particles producing antiprotons and antineutrons have lifetimes just below 1 mm/c. 
Fornengo et al. introduced the condition Ar < 2fm on the spatial separation [11]. In 
Pythia 6.4, this treatment is basically equivalent to considering only antiprotons and 
antineutrons with a production time t = 0 since only non-zero lifetimes that are relevant 
for displaced vertices in collider experiments are stored in the event table. This treatment 
removes all antiprotons and antineutrons coming from metastable mother particles and 
this is the treatment we will adopt for our analysis below. 

not go into the details of the microscopic formation mechanism [29]. In the 1980s it was argued that no 
explicit interaction with a third body is needed since deuterons are formed in a tiny space-time region 
and thus the uncertainty principle applies [43], or that the dominant process for denteron formation 
in a hadronic shower involves a proton and a neutron that are slightly off-shell [44,45]. Nonetheless, 
other authors still stressed the necessity of an explicit interaction with a third body in the nuclear 
shower [46,47]. The authors of [48] recently discussed an alternative denteron formation model that 
makes use of measured cross sections for deuteron production processes with photon or pion emission. 

®We checked how alternative approaches affect the resulting deuteron spectra. Starting from energy 
conservation, i.e. determining the deuteron kinetic energy from the proton and neutron energies and the 
deuteron mass to be T 4 = Ep + — rrid, the resulting spectra differ only notably at deuteron momenta 

below the coalescence momentum. 

^The exact criterion is not stated in the paper; S. Wild, private communication (2014). 


5 




Delahaye&Grefe(2015) 


Delahaye &Crefe(2015) 



Total d yield per Z decay 
— — — onlyAp<po 

Ap <po t = 0 

d yield per had. Z decay in ALEPH range 
onlyAp</7Q 

- Ap </?o “ind t =0 


100 


150 


200 


250 


Coalescence Momentum polMeV) 


Figure 1: Left: Antideuteron yield as a function of the coalescence momentum pQ. The upper 
(blue) lines show the total antideuteron yield and the lower (red) lines show the antideuteron 
yield after applying the same cuts as the ALEPH analysis. The ALEPH result is shown for 
comparison as a grey horizontal band. Right: Antideuteron yield divided by pg as a function of 
Pq. It is clearly visible that the antideuteron yields scale with the third power of the coalescence 
momentum. The solid lines are our main result. The dashed lines show the antideuteron yields 
without excluding antiprotons and antineutrons coming from metastable mother particles. In 
all cases the coloured bands show the uncertainty introduced by Monte Carlo statistics. 


2.2 Determination of the Coalescence Momentum 

Using Pythia 6.4, we simulated 10® events of the process e^e“ —)■ Z to determine po 
from ALEPH data. The authors of [30] found the number of antideuterons produced per 
hadronic Z boson decay in the process —)■ Z —)■ d + X at the Z pole to be® 

Ra = (5.9 ± 1.8 ± 0.5) x 10“® (2.4) 

within the momentum range 0.62 GeV < Pd < 1.03 GeV and in the angular range 
I cosdl < 0.95. In the left panel of hgure 1 we present the simulated antideuteron yield 
as a function of the coalescence momentum. We show the total antideuteron yield per 
Z decay event as well as the yield per hadronic Z decay applying the momentum and 
angular cuts of the ALEPH analysis. The hadronic branching ratio of the Z boson is 
69.91 %. Leptonic Z decays do not contribute at all to antideuteron production [42]. We 
overlay the ALEPH result, adding statistical and systematic errors in quadrature. The 
best-£t values for the coalescence momentum and the l-cx ranges can be basically read off 
this plot and are summarised in table 1 along with the value following from the factorised 
coalescence model, see eq. (2.2). In the right plot of hgure 1 we demonstrate that in all 
cases the dependence of the antideuteron yield on the coalescence momentum is in very 
good agreement with the expected pg behaviour. 

In the left panel of hgure 2, we present the antideuteron spectrum in Z decays as a 
function of the momentum for the three diherent treatments of coalescence we discussed 
above. The spectrum derived from the factorised coalescence approximation is completely 
diherent from the spectra based on event-by-event coalescence. By contrast, the exclusion 

®The first error gives the statistical error, while the second error corresponds to systematic errors. 
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Coalescence Model 

Coalescence Momentum po 

Factorised Coalescence 

141+}^ MeV 

Pythia 6.4 (only Ap < po) 

1731}^ MeV 

Pythia 6.4 (Ap < po and t = 0) 

2031^° MeV 


Table 1: Coalescence momenta for antideuteron production determined from the ALEPH data. 
We compare the factorised coalescence approach using antiproton and antineutron spectra 
generated with Pythia to the event-by-event coalescence using either only the momentum 
condition or additionally excluding antiprotons and antineutrons produced in decays of long- 
lived mother particles. 
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Antideuteron Momentum p (GeV) Gravitino Mass mj/i (GeV) 

Figure 2: Left: Simulated antideuteron spectra from the process e'’'e“ —)• Z —?■ d+X compared 
to the measurement of the ALEPH experiment at LEP. The spectrum obtained from factorised 
coalescence is shown in green. The spectra derived from Pythia are shown in blue (only 
Ap < po) and red (Ap < po and t = 0). The two Monte Carlo spectra are very similar except 
for the high-energy part. The spectrum derived from uncorrelated antiprotons and antineutrons 
exhibits a very different shape. Right: Antideuteron yields in the decay channels Zv, W£ and 
hv as a function of m^i 2 - The yields obtained from the Pythia simulation are independent 
of the gravitino mass and very similar to each other. By contrast, the factorised coalescence 
prescription leads to unphysical behaviour, with yields falling off like 

of antiprotons and antinentrons from metastable mother particles only slightly changes 
the shape of the high-energy part of the spectrnm. The single ALEPH data point cannot 
give any constraints on the shape of the spectrnm. Therefore, all coalescence treatments 
£t the data eqnally well and onr choice in favonr of the Monte Carlo simnlation with 
Ap < Po and t = 0 is only based on the theoretical considerations discnssed above. For the 
rest of the paper we will stick to the event-by-event coalescence nsing po = 203^25 MeV 
and exclnding antiprotons and antinentrons prodnced at displaced vertices. 

Let us conclude this section with a comparison to other po determinations in the 
literature. Our po values differ somewhat from the values found by the authors of [6], 
[10] and [11]. Kadastik et al. used Pythia 8 and found po = 162 ± 17MeV for the 
Ap < Po condition; roughly 10 MeV below our result. Ibarra et al. also used Pythia 8 and 
found a coalescence momentum of po = 192 ± 30 MeV when removing long-lived mother 
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particles; again lOMeV below our result. Fornengo et al. used Pythia 6.4 and found 
Po = 195 ± 22 MeV when excluding metastable mothers. For the Ap < po condition they 
found Po = 180 ± 18 MeV and for the factorised coalescence approximation they found 
Po = 160 ± 19 MeV.® Putze finds po = 202 MeV when using Pythia 6.4, a value similar 
to ours, and po = 195 MeV when using Pythia 8, compatible with Ibarra et al. [49]. 
As mentioned before, differences between different event generators are not completely 
unexpected. However, in some cases there are also different results for the same tool. 
These differences are a bit worrisome since the Monte Carlo statistical uncertainty is 
only of the order of few MeV. A potential source of systematic uncertainty is the choice 
of the Monte Carlo tune, but to our knowledge all of the studies listed above use the 
default tune of Pythia 6.4/8. 

2.3 Antideuteron Spectra from Gravitino Decay 

In models with bilinear i?-parity violation, the gravitino has four main decay channels: 
1 / 3/2 'yu, Zu, W^ and hv [50,51]. As for the case of antiprotons [26], only the channels 
containing a massive gauge or Higgs boson in the hnal state are relevant for antideuterons. 
For the generation of the antideuteron spectra from gravitino decay we simulated 10® 
events with Pythia 6.4 for each of the decay channels and for a set of gravitino masses 
of roughly equal distance on a logarithmic scale: mo /2 = 85GeV, 100 GeV, 150 GeV, 
200 GeV, 300 GeV, 500 GeV, 1 TeV, 2 TeV, 3 TeV, 5 TeV, and 10 TeV. Using the same 
strategy as in [26], we started the Pythia simulation with a resonance decay into two 
particles, Z boson and neutrino, W boson and charged lepton, and Higgs boson and 
neutrino, respectively. In this way the Z, W and Higgs bosons are treated as decaying 
isotropically in their rest frames. 

The antideuteron yields in these decay channels are expected to be independent of the 
gravitino mass since the antideuterons are produced in the fragmentation of an on-shell 
gauge or Higgs boson. Larger gravitino masses should thus only lead to boosted spectra, 
while the underlying physical process remains unchanged.^® Our simulation confirms this 
expectation and shows that the factorised coalescence prescription leads to unphysical 
results, see the right panel of figure 2. The antideuteron yields per gravitino decay are 
summarised in table 2.^^ The resulting spectra are presented in figure 3 for the central 
value of Po = 203 MeV. Although these spectra by eye appear to have the same shape 
as the antiproton spectra from gravitino decay (see [26]), rescaled by a factor of 0(10“'^), 
there is no simple scaling relation among them. 

We will now move on to the discussion of the antideuteron flux observed at Earth. 

® After going back to their code, they now find po = 143 MeV for the latter case, compatible with our 
result; N. Fornengo, private communication (2014). 

^°Note that we did not take into account weak corrections that would lead to a slight increase of the 
yields with increasing gravitino mass, see for instance [52]. 

^^Note that these results differ from those reported in [53] due to an erroneous treatment of the 
Pythia routine in the earlier study. 

^^The decay spectra are available in tabulated form at http://www.desy.de/~mgrefe/files.html. 



Particle type 

Zv 

Wi 

hv 

d + d (t=0) 

(2.15 ±0.69) X 10“^ 

(2.01 ±0.65) X 10“^ 

(3.82 ± 1.22) X 10“^ 

d + d (only p) 

(2.27^i>:S) X 

(2.02^»;i) X 10- 

(4.311};^^) X 10-^ 


Table 2: Multiplicities of deuterons + antideuterons from gravitino decays after the event-by- 
event coalescence process simulated with Pythia 6.4. The central values are given for po = 
203 MeV {t = 0) and po = 173 MeV {p < po), respectively. The quoted errors of roughly 30% 
correspond to the l-u uncertainty in the determination of pq. 
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Figure 3: Deuteron/antideuteron spectra from the two-body decay of a gravitino into Zi' 
(left), Wi {centre) or hv (right). The spectra are shown for the central value of po = 203MeV 
and for gravitino masses of 150GeV (red), 200GeV, 300GeV, 500GeV, 1 TeV (green), 2TeV, 
3TeV, 5TeV, and 10TeV (blue). All spectra are normalized to the respective gravitino mass. 


The flux expected from gravitino decay is simply a linear combination of the fluxes for 
the individual decay channels: 

= BR{Zu) $1^ + BR{Wl) ^ + BR{hp) . (2.5) 

The branching ratios for the different decay channels depend on the choice of the su¬ 
persymmetry parameters and are discussed in detail in [26,53]. Note that they do not 
depend on the amount of i?-parity violation but only on the mass of the gravitino and 
the mass hierarchy of the neutralino sector of the supersymmetric particle spectrum. 
It is hence convenient to label the different cases by the name of the next-to-lightest 
supersymmetric particle (NLSP). In this work we use the values presented in figure 3 
and table 3 of [26]. For illustration, we restrict to an example case where the NLSP is 
Bino-like. 


3 Antideuteron Flux from Gravitino Decay 

As explained in the introduction, antideuterons allow for an almost background-free 
search for an exotic DM component in certain parameter ranges. In fact, no cosmic-ray 
antideuterons have been observed so far and there exists only an upper limit on the 
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antideuteron flux from the BESS experiment [54], In addition, BESS-Polar II looked 
for antideuterons using more than ten times more cosmic-ray data than the previous 
BESS analysis. No candidate antideuterons were observed and a flux limit is expected 
to be published in the near future [55,56]. In addition, several experiments are currently 
taking data or will start operating within the next years to improve this situation: The 
AMS-02 experiment operating on the International Space Station is expected to greatly 
improve on the current sensitivity to antideuteron fluxes [57].^^ Moreover, the General 
AntiParticle Spectrometer (GAPS) is expected to perform several balloon flights, starting 
with a first Antarctic campaign in the austral summer 2019/2020 [21]. A prototype flight 
of the GAPS experiment was successfully carried out in June 2012 [17]. 

Several studies on antideuteron fluxes from DM annihilations or decays can be found 
in the literature [1-5]. These early studies, however, employ the factorised coalescence 
approximation for antideuteron formation, which, as discussed in section 2, is in general 
insufficient to describe the actual production rate. Only more recent studies employ the 
Monte Garlo approach [6-13]. As discussed in section 2, in this work we employ decay 
spectra obtained by this latter method. 

A relativistic antideuteron, formed by the coalescence of an antineutron and an 
antiproton within the hadronic shower of a gravitino DM decay in the Milky Way halo, 
will then propagate through the ISM and might eventually arrive at a detector at Earth. 
As for all cosmic rays, the propagation is described by a diffusion equation of the cosmic- 
ray phase-space density 

V-{V,ij-Ko/3n^ VV') + dE (5ioss ^ - Dee Be^) 

= gp™ + 2 h 5{z) (g^"" + g*""^) - 2 h 6{z) V'- 

For a description of the individual terms we refer to the appendix of [26] and references 
therein. 

3.1 Cross Sections 

The spallation term P®p^^ = nj deserves a separate discussion. This 

term corresponds to the interaction of cosmic-ray antideuterons of velocity with the 
ISM, which is mainly composed of hydrogen and helium. We have considered uh = 
0.9 cm“^ and nne = 0.1 cm“^ [58]. Gontributions of heavier elements have been neglected. 
The annihilating inelastic cross section a™ = is practically given by the 

inelastic cross section since the non-annihilating inelastic cross section is very small due 
to the small binding energy of antideuterons [3,4,9,18]. The inelastic cross section is the 
difference of the total and the elastic cross sections, cr'f®* = — af . 

’ dp dp dp 

A common issue in cosmic-ray antideuteron analyses is that there are no experimental 
data on the cross section of inelastic antideuteron-proton scattering. A typical assump¬ 
tion is hence to rely on the hypothesis that = 2xcrl“^(Tp = Tj/n), where Tj/n 

^^Note that the sensitivity in [57] was estimated based on 5 years of data taking in the superconducting 
magnet set-up of the apparatus. It is not clear whether the actual permanent magnet set-up can reach 
this sensitivity. To date, however, there is no updated study on the AMS-02 antideuteron sensitivity [21]. 
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is the antideuteron kinetic energy per nucleon, see [10] using a direct parametrisation of 
data by Tan and Ng [59] A similar path was followed by [3], calculating the cross 
section as cr^pf{pp) = 2 x {cr^ppipp) — o'ppipp))- Making use of the reasonable assumption of 
symmetry under charge conjugation, one can also use the total antiproton-deuteron cross 
section (J^diPp), for which data exists [42], In this case, the inelastic antiproton-deuteron 
cross section was calculated as crl^®*(pp) = o'^^ipp) — 2 x crpp{pp) [4,9,11]. 

However, although the assumptions presented above give the correct order of mag¬ 
nitude of they are not entirely supported by experimental evidence. In fact, 2 x 
is larger than by roughly 10%. This is theoretically expected due to Glauber 
screening [60] and in nucleon-nucleus collisions one rather expects a geometric scaling 
apA — X app than a linear scaling with nucleon number A [61-63]. 

In this work, we will thus use an inelastic antiproton-deuteron cross section based 
on parametrisations of available data. Indeed, in many previous works it was incorrectly 
stated that there were no data on the inelastic and elastic antiproton-deuteron processes. 
In hgure 4, we compare our parametrisation of with other parametrisations used in 
the literature. We observe that the approach of rescaling by a factor of two over¬ 
shoots the available low-energy data. The parametrisation by Tan and Ng [59] is only 
based on crip®* data from 50 MeV to 100 GeV in antiproton kinetic energy and clearly 
leads to an incorrect high-energy behaviour. The approach of Donato et al. [3], based on 
the difference of crD* and cr|p data that extend to much higher energies, leads to a more 
reasonable high-energy behaviour. Dal et al. [9],^^ using parametrisations of crl^* and crip, 
hnd a better agreement with low-energy data. Our parametrisation gives a comparable 
result, but makes use of better motivated functional forms for the parametrisations and 
explicitly takes into account available low-energy data for the antiproton-deuteron pro¬ 
cess. See appendix A for a detailed derivation of the parametrisation used in this work. 
For crl^“'““ we use the parametrisation of Dal et al. [9]. Our result is certainly not yet 
satisfactory, but we think that it is dehnitely an improvement compared to the treat¬ 
ment of cross sections in earlier antideuteron studies. We hope that our work serves to 
stimulate further discussion in this area. 

In addition to the antideuteron-proton cross section, we need the cross section for 
annihilating antideuteron-Helium scattering. Since no experimental data are available 
for this process, we have rescaled the antideuteron-proton cross sections by the geometric 
factor 4 ^A^ hoping that it suffices to give a reasonable estimate of the true cross section. 

Goncerning the production of tertiaries, for which the differential cross section of the 
non-annihilating inelastic antideuteron-proton process is required, unfortunately no data 
are available. The usual method is thus to take the integrated cross section of the charge- 
conjugate process, and to multiply it by an appropriate energy distribution of the 

outgoing deuteron. In the limiting fragmentation hypothesis, this is simply 1/T', where 

^"^Some earlier studies used the assumption cr|“(rj/n) = 2 x app'^{Tp) [5,7]. This relation, however, 
underestimates since the annihilating antiproton-proton cross section falls quickly with rising 

energy, in contrast to the antideuteron-proton process, see appendix A. 

^^The used parametrisations are not given in the paper; L. Dal, private communication (2014). 
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Figure 4: Comparison of our parametrisation of the inelastic antideuteron-proton cross section 
to other parametrisations used throughout the literature. The data points are taken from the 
Landolt-Bornstein compilation [64], 


T' stands for the kinetic energy of the incoming denteron (see for instance [18,59]). 
However, one can also follow the authors of [3,18] and - inspired by the pp process - use 
the functional form suggested by Anderson et al. [65]: 


dp dX / / . \ n 


non-ann 

pd 


Anderson 

(pO X — -X 

dp 


da 


Anderson 


dp" 


dp" 


(3.2) 


where 


^^Anderson 

dp 



d?a{pp —)■ pX) 
dpdfl 


sin 9 dO 


and 

d?a{pp ^ pX) p^ '-yiE — (3p cos 6) ^ { Vt \ mb 

dp da = ^- E -“P holeel g^' 

In the latter expression, the energy and momentum of the incoming proton are labelled 
E and p, respectively; pt stands for the transverse momentum of the outgoing proton in 
units of GeV. The boost from the laboratory frame to the centre-of-mass frame is ruled 
by the Lorentz coefficients f3 and 7 . We follow this latter method in this work. 


3.2 Propagated fluxes 

As for the case of antiprotons, no primary antideuterons are expected from astrophys- 
ical objects and the dominant background for DM searches are secondary antideuterons 
created in spallation processes of cosmic-ray protons and helium nuclei impinging on the 
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Table 3: Relevant processes for calculating the flux of antiprotons. 
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Table 4: Relevant processes for calculating the flux of antideuterons. 

ISM, i.e. hydrogen and helium gas. In order to estimate the background for the signal 
from gravitino decay, we employ here the same calculation of the astrophysical secondary 
antideuteron flux as used in [3]. As one can see in hgure 5, various processes have to be 
taken into account for the calculation of the antiproton and antideuteron fluxes. The 
relevant processes for antiprotons are presented in table 3 and those for antideuterons 
are presented in table 4. For the antiproton processes, we have used the cross section 
parametrisations given in [66]. For the calculation of secondary antideuterons, we have 
used the antiproton cross sections along with the factorised coalescence prescription for 
antideuteron formation in the hnal state, see section 2}^ 

Since both primaries and secondaries produce tertiaries, i.e. cosmic rays produced 
by the non-annihilating inelastic scattering of high-energy antideuterons on the ISM, we 
have not shown this component separately but rather added it directly to our estimates 
of the primaries and secondaries, respectively. However, we stress that this process, as 
well as convection and diffusive reacceleration, have important impact on the low-energy 
part of the computed fluxes and should not be neglected. Note also that in the case 
of antideuterons the secondaries coming from the spallation of primary antiprotons, i.e. 
antiprotons of DM origin, on the ISM are to be considered as signal and not background. 
This component is much lower than the primary component and the background when 
the gravitino mass is low, but for masses higher than lOTeV, it is of the same order 
of magnitude as the primary component and hence should not be neglected. Figure 5 
displays these different components for a gravitino mass of 100 GeV. 

In order to be able to compare the expected antideuteron fluxes to experimental 
results, one has to take into account the effect of solar modulation [68]. Ideally, one should 
try to fully model this effect, for instance with a simulation like HelioProp [69] as done 
in [11]. However, this requires knowing precisely the status of the solar environment at 

Alternatively, one could also use the Monte Carlo approach for calculating secondaries [67]. 
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Top of the Atmosphere Fluxes 



Figure 5: Various components of the antiproton and antideuteron fluxes at Earth. The as- 
trophysical background (secondaries) comes from the interaction of cosmic-ray protons and 
a-particles with the ISM. For the case of antideuterons also interactions of secondary anti¬ 
protons with the ISM contribute. The primary antiproton and antideuteron components come 
from gravitino decay; for the case of antideuterons also from the interaction of primary anti¬ 
protons with the ISM. The tertiary component corresponds to a redistribution of high-energy 
cosmic rays to lower energies due to non-annihilating inelastic scattering off the ISM. Since 
every component creates tertiaries, they have been incorporated directly and are not displayed 
as separate components. 

the time of data taking. In the absence of data and because we do not know yet the time 
of the data taking, we think it is untimely to go through such a precise modelling. For 
Figure 5, we hence satisfy ourselves with the so-called Fisk approximation [70]: 

4|0A(£T0a, ^ ^ gXOA ^ _ (3 3, 

assuming a Fisk potential of (pp = 500 MV as an example since this values allows to have 
a good agreement with the most recent PAMELA antiproton data [71]. In the subsequent 
hgures, where we do not show antiproton data, we will only consider interstellar fluxes 
since we do not know what the Fisk potential will be when data will Anally be taken. 

In flgure 6, we present the interstellar antideuteron spectrum from gravitino DM 
decays and compare it to the expected astrophysical background and the flux limit 
obtained by the BESS experiment [54]. In addition, we present the projected sensitivity 
regions of the BESS-Polar II [56], AMS-02 [57] and GAPS [17]^”^ experiments. In the 
left panel, we show the uncertainty band due to the lack of precise knowledge of the 
propagation parameters as constrained by measurements of the boron-to-carbon ratio in 

^^The GAPS sensitivity assumed in this work corresponds to three Antarctic long duration balloon 
flights with a total duration of 105 days [21]. 
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Propagation Uncertainties 


Antiproton Constraints 



Figure 6: Left: Cosmic-ray antideuteron flux expected from the decay of gravitino DM com¬ 
pared to the expectation from astrophysical secondary production and the sensitivities of forth¬ 
coming experiments. The flux from gravitino decay is shown for a lifetime of 10^® s and masses 
of 100 GeV, 1 TeV and lOTeV. The coloured bands correspond to propagation uncertainties 
within constraints from boron-to-carbon ratio measurements. Right: Same as left panel but 
fixing the gravitino lifetime to the lowest value allowed by antiproton constraints. No solar 
modulation has been implemented here since the time of data taking is unknown. 


Model 

5 

Ko (kpc^Myr) 

L (kpc) 

Vc (km/s) 

MIN 

0.85 

0.0016 

1 

13.5 

MED 

0.70 

0.0112 

4 

12 

MAX 

0.46 

0.0765 

15 

5 


Table 5: Parameters of cosmic-ray propagation models that correspond, respectively, to the 
best fit of cosmic-ray boron-to-carbon data (MED) as well as the minimal (MIN) or max¬ 
imal (MAX) antiproton flux compatible with cosmic-ray boron-to-carbon data. Figures taken 
from [73]. 


cosmic rays [72] . The coloured bands correspond to a full scan over the allowed parameter 
space. As an illustration we also display the fluxes obtained with three benchmark models 
often used in the literature (see table 5). We use a common decay lifetime of T 3/2 = 
10^® s for illustration. Note that, as for the case of antiprotons, the MIN/MED/MAX 
benchmark models do not size the full extent of the uncertainty band. This shows again 
the importance of performing scans over the full propagation parameter space allowed 
by other data rather than checking only a few cases. 

In the right panel of figure 6, we display the same fluxes but setting the decay 
lifetime to the minimum values allowed by the constraints obtained using the antiproton 
measurements (see [26]). Note that the coloured bands correspond to the extreme cases 
obtained for antiprotons and are not the full uncertainty band. Indeed, the limiting cases 
correspond to situations where the antiproton flux is maximised at a given energy bin. 
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Antiproton Constraints 



Figure 7: Same as the right panel of Fig. 6 for the MED propagation parameters. The only 
difference is that the fluxes displayed are now corrected for solar modulation (cpp = 300 MV 
and 4>f = 700 MV). As one can see, because the fluxes are relatively flat below lOGeV, solar 
modulation has little impact on the expected fluxes. 


This does not mean that the antideuteron flux is maximal over the whole energy range. 

As one can see from Fig. 7, if the Fisk approximation describes correctly the influ¬ 
ence of solar modulation on antideuterons, then this should not affect our conclusions 
dramatically. Indeed, unlike many other cosmic-ray species, antideuteron fluxes are rel¬ 
atively flat below 10 GeV and a shift of the antideuteron energy does not affect the flux 
very strongly, at least within a reasonable range of the Fisk potential. 


4 Discussion of the Detection Prospects 

A striking feature of figure 6 is that for masses as low as 100 GeV the gravitino decay 
signal can be of the same order as the astrophysical background below a few GeV, even 
for lifetimes as large as 10 ^® s, a value not yet excluded by gamma-ray and antiproton 
observations (see for instance [74]).^® In this respect, it would also be interesting to see 
what antideuteron flux could be expected for even lower gravitino masses. It could thus 
be worthwhile to study this region in a future work using the spectra obtained from 
gravitino three-body decays [53]. 

But also for larger gravitino masses the antideuteron signal could be at the same 
order as the background. This was not observed in earlier studies as the signal for large 
DM masses is artificially suppressed in the factorised coalescence prescription. Therefore, 

^®See also [75, 76] for recent works on gamma-ray constraints on the gravitino lifetime, taking into 
account the latest Fermi LAT data [77]. 
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for decaying DM candidates there is in principle also the possibility of observing an exotic 
component in the higher-energetic part of the spectrum, where currently no experiments 
are planned. Note, however, that both background and signal are extremely low and 
quite challenging for experimentalists as this would mean improving sensitivity by at 
least four orders of magnitude in flux, but also to reach much higher energies. 

When taking into account the constraints derived from antiproton observations [26], 
we hnd that the remaining parameter space for having a gravitino decay signal signihc- 
antly higher than the astrophysical background becomes quite small but does not vanish 
completely. Since the coalescence process still suffers large theoretical uncertainties, one 
cannot exclude that all the fluxes are in fact larger (or smaller) than what we assume 
here. The l-cx uncertainty in pq could lead to an increase of roughly 30% in the antideu- 
teron flux from gravitino decay, irrespective of the antiproton constraints. In addition, 
making use of Monte Carlo methods to estimate the production of secondary antideu- 
terons instead of the factorised coalescence model used here, tends to predict a slightly 
lower background level at low energies [67]. This would increase the signal-to-background 
ratio, independently of any other constraints. 

Still, it seems clear that the current generation of experiments will not be able to 
observe any antideuteron events. Only an improvement of the flux sensitivity by two 
orders of magnitude should at least allow for a detection of astrophysical antideuterons 
- or even those coming from gravitino decay. The main hopes seem to reside either in 
the highest energy range (above ~ 50 GeV) or in the lowest one (below ~ 1 GeV). Note, 
however, that the latter is affected by solar modulation, a phenomenon that to date is 
not fully under modelling control. 

This clearly challenges antideuterons as the golden channel it has long thought to 
be. Indeed the antiproton channel has become extremely constraining thanks to the 
PAMELA data [71]. A forthcoming release of AMS-02 antiproton data could make these 
constraints a bit more stringent, especially at higher energies. 


5 Conclusion 

In this work, we have studied the potential of detecting cosmic-ray antideuterons pro¬ 
duced in the decay of gravitino dark matter within a framework of bilinear i?-parity 
violation. This work was a natural sequel of a related work concerning cosmic-ray an¬ 
tiprotons. After discussing the deuteron formation in hadronic showers and calculating 
antideuteron spectra from gravitino decay, we have determined the gravitino decay sig¬ 
nal at Earth. We have also assessed the uncertainties affecting the expectations for the 
astrophysical background and the signal. We have shown that there is some room left for 
a discovery of gravitino decays through antideuterons, however not within the sensitivity 
of the current and planned generation of experiments. 

We have also shown that - once the antiproton constraints are taken into account 
- the remaining parameter space for a detection of gravitino dark matter with bilinear 
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i?-parity violation via antideuterons is quite small. On the other hand, since not much 
progress is expected in the background-limited antiproton channel in the coming years, 
the antideuteron channel could still serve to put stronger constraints on the strength 
of the i?-parity violation in the future. If the detection technology were to improve 
considerably, also the high-energy regime (above ~ 50 GeV) would become interesting 
for the search of gravitino dark matter. 
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A Cross Sections 

In this appendix, we briefly present the parametrisations for the antiproton-proton and 
antiproton-deuteron cross sections we used to estimate the inelastic scattering cross 
sections relevant for the spallation term of antideuteron propagation in the Milky Way. 


Antiproton—Proton Cross Sections For the total cross section of antiproton-proton 
scattering, a large number of data points ranging from roughly 200 MeV to 2PeV in 
antiproton momentum (in the rest frame of the proton target) exist in the literature [42]. 
A useful parametrisation is given by [78,79]: 




asmpt V 


1 + 


^s-4m^R|](s) 



di I ^2 da \ 


(A.l) 
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where s = 2 rrip + 2 mpyrrip + Pp is the centre-of-mass energy of the pp system and 
'^asmpt('5) = 00 + ^2 ln^(\/s/\/^) 5 a/^ = (20.74 ± 1 . 21 ) GeV, 

(jtot ('e') 

Rl{s) = 2.568 GeV^ mb“^ ^ - B{s ), B{s) = &o + &2 

ZTT 

In the dehnition of Rq, the factor Imb = 2.568 GeV~^ simply accounts for unit con¬ 
version. The remaining parameters in these expressions were determined from a £t to 
experimental data in [78]: 

ao = (42.05 ± 0.11) mb, aa = (1.755 ± 0.083) mb, 

6 o = (11.92 ±0.15) GeV-^ 62 = (0.304 ± 0.019) GeV-^ c = (6.7 ± 1.8) GeV^ 

di = (-12.1 ±1.0) GeV, da = (90 ± 16) GeV^ 4 = (-111 ± 22) GeVl 


With these parameters, eq. (A.l) gives a relatively good fit to the available data. Includ¬ 
ing also high-energy proton-proton data between 10 TeV and 2 EeV in proton momentum 
- where the antiproton-proton and proton-proton cross sections should be equivalent - 
we get a goodness of fit of y^/dof = 2821/460.^® Given that many early experimental 
works, especially from the 1960s and 1970s, only quote statistical errors without any 
assessment of systematic uncertainties, the £t is quite acceptable (see also the remark 
on errors in the introduction of [64]). 

For the elastic cross section of antiproton-proton scattering, there is also a lot of data 
ranging from roughly 200 MeV to 2 PeV in antiproton momentum [42]. Using typical basis 
functions for cross sections (see [64]), we fit the data using the piecewise ansatz:^° 



Oq ± ai(pp/GeV)" if Pp < 5 GeV 

bo ± bi{pp/GeV)^ ± 62 ln(pj5/GeV) ± 63 ln^(pp/GeV) if pj, > 5 GeV 


and requiring differentiability at Pp = 5 GeV. We find the parameters 


ao = —35.71mb, ai = 81.03mb, n =—0.2567, 60 = 10.49 mb, 

5i = 86.04mb, m =—1.385, 62 =—1-360 mb, 63 = 0.1312mb, 

giving a goodness of fit of y^/dof = 364.6/167.^^ Besides data on elastic antiproton- 
proton scattering we included data on elastic proton-proton scattering with incoming 
momentum between 100 GeV and 34 PeV [42]. 

For the inelastic antiproton-proton cross section there is considerably less data than 
in the previous cases, ranging only from 300 MeV to 175 GeV [64]. This cross section is, 

have also checked that the total cross section data for the antiproton-neutron [42] and 
antineutron-proton [64] scattering processes are compatible with the antiproton-proton result. 

^°The authors of [80] present a parametrisation of the elastic antiproton-proton cross section based 
on the methods of [78,79]. However, it appears that their parameters are faulty since we were not able 
to reproduce their result. 

^^We have also checked that the cross section data for the elastic antiproton-neutron [42] and 
antineutron-proton [64] scattering processes are compatible with the antiproton-proton result. 


19 





Delaha^e & Grefe (20iS) Delahaye & Grefe (2015) 



Figure 8: Left: Total, elastic, inelastic and annihilating inelastic cross sections of the pp 
process. Data points are taken from the Particle Data Group [42] and the Landolt-Bornstein 
compilation [64]. The lines are the parametrisations discussed in the text. Right: Same as left 
panel but for the cross sections of the pd process. 

by definition, the difference between the total and the elastic cross section and we just 
check the consistency with available data. Including also data on inelastic proton-proton 
scattering between 1 and 2 TeV, we hnd a goodness of £t of x^/dof = 404/43. Given that 
only three of the 43 data points come with an estimate of systematic uncertainties, we 
think this is an acceptable result. 

An overview of the antiproton-proton cross section data is presented in hgure 8 . In 
addition to the total, elastic and inelastic cross sections, we present data on annihilating 
inelastic antiproton-proton scattering ranging from 240 MeV to 22 GeV [64] . These data 
can be parametrised by the function 

(Tpp^iPp) = 0.2367mb + 63.86 mb (pp/GeV)“°'®®®, (A.3) 

giving a goodness of fit of y^/dof = 46.5/27. 

Antiproton—Deuteron Cross Sections The total cross section for antiproton- 
deuteron scattering is almost twice as large as the total antiproton-proton cross sec¬ 
tion. The exact result depends on the size of the elastic and inelastic Glauber shadow or 
screening corrections [60,81,82]: 

^tot o ^tot rel jrinel 
(^pd =^(^pp - ^ ■ 

If these correction terms are known, also the elastic and inelastic cross section components 
can be easily determined individually: 

^el o rel ^inel o rinel 

(^pd = ^(^pp-o , ap^ =2 app -d . 

Arkhipov describes antiproton-deuteron scattering and the corresponding correction 
terms in [82]. Unfortunately, his parametrisations do not describe very well the avail¬ 
able data. Therefore, we will follow an alternative route. 

have also checked that the inelastic antiproton-neutron scattering data [64] are compatible with 
the antiproton-proton result. 
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For the total cross section of antiproton-deuteron scattering, there is only data from 
roughly 300 MeV to 280 GeV in antiproton momentum [42] . For the high-energy part we 
thus assume that the cross section matches the antiproton-proton cross section, rescaled 
with a suitable factor. We find that eq. (A.l), rescaled with a factor of 1.85, gives a 
reasonably good fit to the antiproton-deuteron data above 6 GeV. Using the same basis 
functions as for the elastic antiproton-proton cross section, we then fit the data using 
the piecewise ansatz: 


^pd (Pp) = 


ao + ai{pp/GeVy 


if Pp < 1-35 GeV 


^0 + ^i(Pp/GeV)’" -|- 62 lii(pp/GeV) -I- 63 ln^(pp/GeV) if pp > 1.35 GeV 


(A.4) 

and requiring differentiability at pp = 1.35 GeV. Including rescaled antiproton-proton 
data above ITeV and rescaled proton-proton data above lOTeV in the £t, we hnd 


00 = 75.59 mb, oi = 134.3 mb, n =—0.6647, 60 = 77.33 mb, 

61 = 132.4 mb, m =—0.6388 , 62 = —3.638 mb, 63 = 0.5615 mb, 

giving a goodness of £t of y^/dof = 452.9/159. 

For inelastic antiproton-deuteron scattering, there are data from roughly 300 MeV 
to 2 GeV [64]. As for the case of antiproton-proton scattering, we will determine the 
inelastic cross section by subtracting the elastic from the total cross section. Lacking a 
physically motivated ansatz for the elastic antiproton-deuteron cross section, we rescale 
eq. (A.2) by a suitable factor to match the data points. We find that 

^pfiPp) = ^IdiPp) - 1-75 X afpipp) (A.5) 

gives a goodness of fit of y^/dof = 15.8/13. We still have to check if our approach for 
the elastic cross section matches the available data ranging from roughly 300 MeV to 
5 GeV [64]. As can be seen from the right panel of hgure 8 , the shape of the rescaled 
elastic cross section does not entirely follow the data points. The goodness of fit thus is 
a very poor y^/dof = 535/9.^^ Gonsidering only the low-energy data would result in an 
acceptable y^/dof = 46/9. 

For the non-annihilating inelastic antiproton-deuteron cross section that is necessary 
for the calculation of tertiaries we use the parametrisation of [9]:^'^ 

^non-anin'p_^ _ j^g-2.141+5.865exp(-logi|3(pp/GeV))-3.398exp(-2logiQ(pp/GeV)) (A 6 ) 
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